Many-body correlations in Semiclassical Molecular Dynamics and Skyrme interaction 
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Constraint Molecular dynamics CoMD calculations have been performed for asymmetric nuclear 
matter (NM) by using a simple effective interactions of the Skyrme type. The set of parameter 
values reproducing common accepted saturation properties of nuclear matter have been obtained 
for different degree of stiffness characterizing the iso-vectorial potential density dependence. A 
comparison with results obtained in the limit of the Semi-Classical Mean Field approximation using 
the same kind of interaction put in evidence the role played by the many-body correlations in 
to explain the noticeable differences obtained in the parameter values in the two cases. Even 
if from a numerical point of view the obtained results are strictly valid for the CoMD model, 
some rather general feature of the discussed correlations can give a wider meaning to the obtained 
differences being strongly related to the spacial correlations generated in the semiclassical wave 
packets dynamics. 



PACS numbers: 24,10.-i,24.10.Cn, 21.65.Ef, 21.65.Mn 



I. INTRODUCTION 

The description of many-body systems is one of the 
most difhcult problems in nuclear physics due to the 
complexity of this kind of systems which are quantum 
objects described by a large number of degrees of free- 
dom. A large variety of theoretical models have been 
developed using mean-field based and beyond-mean field 
approaches like Density Functional Theory I'l and En- 
ergy Density Function Theory j2|, i3|] . In these approaches 
which use the independent particle approximation as a 
starting point, phenomenological effective interaction like 
Skyrme and Gogny forces are widely used taking ad- 
vantage of their simple formQ. In nuclear structure 
modeling we can quote the Skyrme-Harthree-Fock (SHF) 
methods, the relativistic mean-field (RMF) approach (as 
well as their Bogoliubov extensions), and the Harthree- 
Fock-Bogoliubov methods with the finite-range Gogny 
force 0. As an example we can quote the work of P.- 
G.Reinhard et al fo"] in which the liquid drop parameters 
have been obtained starting from microscopic SHF and 
RMF calculations with a rather wide variety of Skyrme 
interactions. Each of them produce energy functionals 
with parameters adjusted to reproduce the basic nuclear 
bulk properties in the valley of stability. 

Complexity still becomes higher when nuclear dynam- 
ics, triggered by nuclear collisions, is studied to under- 
stand the property of nuclear forces far from the stability. 
At the Fermi energy and beyond, semiclassical methods 
become necessary to describe the produced processes in 
which practically all the degrees of freedom are involved. 
Many-body correlations are responsible for the reorga- 
nization of the hot Nuclear Matter(NM) in to clusters. 
Two main classes of approaches have been developed to 
handle this complex scenario. The first one is based on 
the Boltzman transport equation including two and, at 
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the higher energy, also three body collisions term0-Q. 
In this case the theoretical approach is able to describe 
the time evolution of the one-body distribution function 
in phase-space. Apart from the collision term, the Boltz- 
man equation corresponds to the semi-classical limit of 
time-dependent Harthree-Fock equations in phase-space. 
Different models have been implemented on this ground. 
They essentially differ from each other in the strategy 
adopted to simulate the collision terms and in the way 
to represent the phase-space distribution starting from 
the nucleonic degree of freedom (test particle methods). 
As in the Harthrre-Fock mean field method the main in- 
gredient concerning the interaction is the energy density 
which for a Skyrme type phenomenological interaction 
has a rather manageable form. In particular the local 
part of the interaction produces a simple functional for 
the potential energy expressed through an algebraic func- 
tion of the density (see also Sec. II). In the early 90 a 
further development included stochastic forces (typical 
of Langevin processes) in the so called Stochastic 
Mean Field model to produce fluctuations capable 
of describing associated phenomena like mean-field in- 
stabilities leading to the cluster production. The sec- 
ond line of development of theoretical semiclassical ap- 
proaches to deal the nuclear many-body problem is rep- 
resented by the so called Molecular Dynamics models. 
The point of view of these methods is some sense an- 
tithetic compared to the ones based on the mean-field 
concept. However also in these cases, due to their sim- 
plicity the same kind of phenomenological effective inter- 
actions Skyrme and Gogny are widely used. The start- 
ing point of all these models is a fundamental assump- 
tion on the wave function describing the single nucleon 
which is represented through a wave packet well local- 
ized in phase-space with uncertainty satisfying the un- 
certainty principle. The many-body wave fucntion can 
be expressed through a direct product, leading to the so 
called Quantum Molecular Dynamics (QMD)-like mod- 
els [13 - [l5| or an anti-symmetrized product like in FMD 
and AMD [H \vf\ . Variational principles have been used 



2 



to obtain the equation of motion for the many-body 
wave functions. Within their own approximation schemes 
these models treat the many-body problem without the 
usage of the mean-field approximation. Many-body cor- 
relations are spontaneously produced and are able to de- 
scribe the main features concerning the multi-breakup 
processes [l^ observed in heavy ion collision at interme- 
diate energy and beyond. Due to these deep conceptual 
differences between the mean-field based models and the 
molecular dynamics ones it becomes interesting, and nec- 
essary in our opinion, to investigate the differences be- 
tween the energy-density functionals which are produced 
with these two classes of approaches when one uses, in 
both the cases, the same kind of effective elementary in- 
teraction between nucleons. The study presented in this 
work is performed by us ing the Constraint Molecular Dy- 
namical Model CoMD [Tsl . [T^ , and the comparison will 
concern both the saturation properties of the symmetric 
NM and the properties of the symmetry energy obtained 
with a simple Skyrme force. The study is limited to mod- 
erate values of the asymmetry parameter as the ones re- 
ally encountered in nucleus-nucleus collisions. As it will 
be shown the presence of iso-vectorial forces strongly af- 
fects the results of this comparison (see also [111). The 
work is organized as follows: in Sec. II we illustrate the 
choice of the effective interaction and the related total en- 
ergy functional are introduced. In Sec. Ill we describe the 
NM simulations. The results are presented in Sec. IV. 
Finally Sec.V contains the summary and the concluding 
remarks. 



II. THE EFFECTIVE INTERACTION AND 
THE TOTAL ENERGY 

Before illustrating with some detail the model and the 
NM simulations, in this section we briefly introduce the 
effective interaction which we use in our calculations. 
We also present the way in which the related total en- 
ergy is obtained in the case of the Semiclassical Mean- 
Field approximation (Se-MFA) and in the case of the 
semiclassical- wave packet dynamics. 

The elementary interaction between two nucleons with 
spacial coordinates r,r' and third isospin component t,t' 
is of the Skyrme type and has the following form: 

F(r,r') = V^^^H{v-y') = 
To 



{a + l)p°. 



S{r - r') 



(1) 



+ —Fl^{2Srr'-l)Sir-r') 
Pa 

The first term of cq(l) represents the iso-scalar contri- 
bution to the two-body interaction, the second one is 
the usual 3-body effective interaction. For this term the 



p density dependence is modeled through the parame- 
ter a. This term appears to be a generalization of the 
Brink- Voutherin Q three-body term corresponding to 
a — 2. This generalization is sustained both on the ba- 
sis of phenomenological parametrization used in mean- 
field approaches [11 and on more complex microscopic 
calculations based on G matrix calculations with realis- 
tic interactions. The third term represents instead the 
Iso- Vectorial contribution. The form factors Fk have the 
following expression: 
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Po is the saturation density value. Also in this case the 
F form factors are introduced to reproduce the results 
of more complex microscopic calculations based on real- 
istic interaction concerning the symmetry energy and re- 
flecting effects beyond the two-body interaction (20l - [23 |. 
These form factors have been widely used in Hartree- 
Fock calculations, in semiclassical BUU calculations and 
in several molecular dynamcs models. In particular, F4 
will be used in the limit p/po ^ 1 to perform calculations 
with stiffness parameter 7 different from the ones related 
to the others form factors (in the same limit Fi , F2 , F3 
correspond to 7 values 1.5, 1 and 0.5 respectively). Fi- 
nally we observe that for simplicity reason we do not add 
to this interaction non local effects. Regarding this as- 
sumption we also note that some of the suggested Skyrme 
parametrizations have a vanishing or very small terms 
related to the non locality at the saturation density in 
the limit of Nuclear Matter (see for example Ref. 25]). 
Moreover, as it will be shown in the following, in CoMD 
calculations the effective potential which determines the 
effective force able to govern the motion of the wave pack- 
ets centroid is already non-local (sum of Gaussian con- 
tributions depending on the intra nucleons distances) . It 
in fact represents the result of the convolution of the well 
localized wave packets with 6 function in space. 

Even if well-known, we shortly recall in the next section 
the main ingredients of the Se-MFA. We think this sec- 
tion useful to present in the following the way in which 
this assumption is instead broken in the semiclassical 
wave-packets dynamics. 



II. 1. The Semiclassical Mean Field Approximation 

At a flxed time starting from eq.(l) (we are consid- 
ering here a stationary problem) , in a rather general 
way, we can obtain the expression for the total energy W 
related to the potential interaction by folding the elemen- 
tary interaction with the two-body density distribution 
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in phase-space D{r, r', p, p'). By taking into account the 
usual truncation condition on D typical of the mean-field 
approximation, we can write: 



D = D,ir,p) ■ D,ir\p') 



(7) 



where Di is the one-body distribution. Therefore by us- 
ing eq.(l) for W we get: 



W^^J V{r,r')D{r,r',p,p')drdr'dpdp' 
= V^^^D,{r,p)Di{r,p')drdpdp' 



(8) 



Taking in to account that by definition J Di(r,p)dp — 
p{r) we obtain: 



(9) 



i]/(2)p2 appear to be the energy density associated to 
the potential energy from which we can obtain the re- 
lated binding energy per nucleon due to the potential 
interaction Epot 



pot — Utwb + Utrb + Uasy 
P 



2 2 po (a- + 2 



(10) 



/3 



represents the charge/mass asymmetry pa- 



rameter evaluated from the neutron and proton density 
Pn , Pp respectively. The total binding energy E is ob- 
tained by adding the kinetic contribution coming from 
the Fermi motion 



Resonances [26l - l28j . Finally we note that the structure of 
the obtained energy- density functional in Se-MFA makes 
independent the choice of the parameter values describ- 
ing the symmetry energy from the other ones which in- 
stead are fixed from the saturation properties of the sym- 
metric nuclear matter. We will show in the next section 
that this is no longer true in the CoMD approach. 



II. 2. The case of the CoMD model 

In molecular dynamics approaches a basic assumption 
is done on the wave functions describing the nucleonic de- 
gree of freedom. It is commonly assumed that the wave 
function is a gaussian wave packet with fixed width (7^ 
in coordinate space. The centroid in phase-space is indi- 
cated with Vi, Pi 



-exp[ 



(r 



.rpi 



(2^(72). 

The Wigner transform of is 



2ai 



{2Trara-py 



7exp[- 



(P- 



2cr2 



2al 



(15) 



(16) 



the widths in momentum and space satisfy the minimum 
uncertainty principle condition (T^Cp = 5^. Another 
assumption concerns the N-body Wigner distribution 
which in CoMD model ^](like in QMD approach[12]) 
is a direct product of the single particle distributions. In 
this case therefore for a system formed by A nucleons 
the one-body and 2-body distributions above introduced 
have a multi-component structure that is: 



E - Epot + £^fi„ (11) 

In , (3 terms of order greater than two are neglected. 
In particular we see how the iso- vectorial vectorial forces 
with strength proportional to T4 contribute to the sym- 
metry energy Easy depending on the asymmetry param- 
eter p. For quadratic form in /3 we get: 



l,d'^E 
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(13) 



(14) 



i^t^Hp) 



The common accepted bulk value of esym at the satura- 
tion density is about 30 MeV even if relativistic Hartree- 
Fock models can predict higher values up to about 40 
MeV [2^ [13. However, still under study is the den- 
sity dependence of this quantity which is able to af- 
fect neutron-skin thickness in nuclei and Giant Monopole 



Di(r,p) =^/,(r,p) 



(17) 



i?(r,p,r',p')= /«(r,P)/.(r',p') (18) 

From the above relations we can see that in general 
D ^ _Di(r, p)L)i (r', p'). In particular, for the case in 
which these distributions are expressed as a sum of differ- 
ent localized components inside a volume Vg , we can in- 
dicate with Qg the number of components which give non 
negligible contributions in Vg. The relative difference as- 
sociated to eq.(18) 1 — Di(r,p)Di(r ,p ) estimated to 
be of the order of 1 /ag which corresponds to the ratio be- 
tween diagonal (i = j) and off-diagonal elements {i ^ j) 
within the ensemble of ag{ag — 1) terms. In semiclassical 
mean-field models ag can be enough high to make the 
difference negligible, in fact the single particle distribu- 
tion usually spreads over the whole system (test particles 
methods) and therefore the truncation condition eq.(18) 
can be retained valid 0-0 • On the contrary this is not 
surely the case for the molecular dynamics approaches for 
which the typical spreading volume if of the order of 2-10 
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fm^. Localization and therefore coherence of the wave- 
packed used to describe the single-particle wave- functions 
allows to keep memory of the two-body nature of the 
inter-particles interaction, and at same time, allows for 
the spontaneous appearance of the clustering phenomena 
in simulations concerning low-density excited portion of 
nuclear matter. With these assumptions on the 2-body 
phase-space distribution, taking in to account the proper- 
ties of the 6 function, we can obtain the explicity expres- 
sion for the different terms concerning the total energy; 
for the two-body isoscalar contribution Wtwb we get: 



we obtain: 



Wtwb = 



Wtwb = 
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(21) 



In the above expression SI is the normalized sum of 
the Gaussian terms and it represents just a measure of 
the overlap between the nucleonic wave-packets. Its two 
body character is quite explicit. In the calulations con- 
cerning the NM simulation that we will illustrate in the 
next sections , the large number of particles A involved 
in the system allow us to write the above quantity in a 
simpler way by introdcing the average overlap per nu- 
cleon Si = Sy whose dependence on the particle index in 
the ideal case can be omitted: 



Wt,„h = 



twb 



Wtwb 



2po 



Sy 



(22) 
(23) 



By comparing the expression obtained for Etwb in the two 
different appraches (eq.(lO) and eq.(23)) we note a formal 
analogy where the variable p is substituted by the over- 
lapp integral Sy per nucleon. We however observe that 
this analogy is only formal, this aspect will be discussed 
in some detail in the next subsection. Concerning the 
three-body term, according to the evaluations reported 
in reference (T2| and taking in to account the previous 
obervations we get: 



Etrb — 
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(a + 1)pE^ 



(24) 



For the term related to the iso-vectorial interaction and 
for the most simple case F' = 1 in the limit A^N^Z » 1 
(N and Z represents the number of neutron and protons 
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(25) 
(26) 

(27) 
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where p'^'^ with cc' equal to nn, pp and np represents 
the overlap integral per couples of neutrons, protons and 
neutron-proton. A more convenient form for the above 
expression is obtained by introducing the two followng 
quantities: 



N^ + Z^ 
'ff^-p 



taking into account that by definition 
Z = - — we get: 



W^ - — 
2po 

E?s. = ^/k(S:)pA[il + ^)0 

Ebias = --. — Fk{Sv)pAa 
4po 



(29) 
(30) 



and 



(31) 
(32) 
(33) 



-] 
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where pA = Ap. The expressions in eqs. (32,33) also con- 
tain a generalizzation to the cases in which we use the 
generic form factors F^. Here i^^ keep the same func- 



tional form in 5*1, using the formal analogy ^ 



above discussed. 

From eq.(32) we obtain that for o; ^ (as CoMD 
calculations predicts, see next sections) the iso-vectorial 
force produces, beyond a modified term for the symme- 
try energy, also another term (the second term) which 
we name -Ehias (eq.(33)). The effect of this last term is 
not-negligible if compared to the balance of the different 
term appearing in the expression of the total energy. The 
correlation coefficient a by definition (see eq.(30)) repre- 
sents the difference in percentage of the overlap between 
the neutron-proton couples from the one related to the 
couples formed by homonym nucleons. It also depends 
on the strength of the iso-vectorial forces [Tj^ parameter) . 

Finally, according to the general strategy character- 
izing the CoMD approach [H, [13], the kinetic energy 
contribution E^in is obtained in a self-consistent way 
through the numerical constraint. Minimum energy con- 
figurations are obtained by applying the cooling- warming 
procedure coupled with the constraint on the occupation 
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FIG. 1. The average kinetic energy per nucleons Ckin is shown 
as a function of the reduced density for a symmetric system 
containing 3500 particles. In the panel (a) we show the case of 
the non-interacting case and in the panel (b) a typical case in 
presence of the inter-particles interaction. The continue lines 
represent the p^^^ trend as expected for Fermionic systems 
(see the text). 



numbers related to the the single particle wave functions 
as required by the Pauli principle. As an example in 
Fig.l we show ekm = for a symmetric systems con- 
taining 3500 particles at different densities. The results 
have been corrected for surface effects (see Sec. III. 1). In 
panels a) and b) we show with black points the case of 
non-interacting and interacting particles respectively. In 
both the cases, within the uncertainty of the numeri- 
cal procedure, the behavior as a function of the density 
is well represented by a p^^^ dependence typical for the 
fermionic systems. This is shown in the figure trough 
slight continuous lines. In particular, the above men- 
tioned density dependence well describe also the case of 
an interacting system of particles. In fact, as previously 
discussed, the average field regulating the related dynam- 
ics is defined (see eqs. (22-33) over a rather large number 
of particles and it is performed after having included the 
two-body particle interaction. In the following (sec sec- 
tion III.l) we will represents the kinetic energy contribu- 
tion through a p'^^^ density dependence best-fitting the 
model calculations. The total energy per nuclcon is there- 
fore obtained by adding all the contributions discussed in 



this section and it will be indicated as = Epot{Sv, a, 
Pa,P) + Ekin{p,l3). In particular we note that, at this 
level, E'-' depends on the new defined primary quantity 
Sv, a, pA 

In the following section we will try to relate these quan- 
tities to more fundamental ones that are the density p and 
the spatial correlation function v between nucleon pairs. 



II. 3. Overlap integrals and spacial correlations 

In a rather general way and as suggested from CoMD 
calculations (see Fig. 2), for an uniform many-body sys- 
tem at a fixed density we can introduce the probability p 
to have a particle in the volume dVi localized in ri and 
a second one in the volume dV2 localized in r2. Due to 
the uniformity condition p depend only on r = |ri — r2|. 
p can be expressed in the following form: 



p=l± kov{r) !/(0) = 1 fco > 



(34) 



Moreover \mir^^v{r) = 0, that is: no spacial correla- 
tions can be expected for very distant particles. Non 
zero values of v can be instead expected at relatively 
small distances due both to the interaction (for an at- 
tractive interaction the positive sign must be considered 
in eq.(35)) and to the symmetry of the many-body wave 
function for quanta! systems of identical particles (in the 
case of identical Fermions we must consider the sign mi- 
nus and fco = 1 )• For a classical system of non interact- 
ing particles we have a vanishing correlation effect p=l 
(fco = 0). 

For our aims, we need to evaluate a normalized prob- 
ability P = cpin such a way: 



/ PdV = 4ttc [ 
Jv Jo 



p{r)r^dr — A 



V-Vc 



(35) 



Vc represents the volume in which the well localized cor- 
relation function v is different from zero. 

This volume is always finite and of the order of the <Jr^ 
in our model calculations. So that in the limit V ^ oo 
we get c = p and : 



P(r) = p(l ± fcoz^(r)) 



(36) 



In Fig. 2, just as an example, we display results of cal- 
culations for a given set of parameters for the Skyrme 
interaction (see the next section) . The calculations show 
the probability to find two nucleons at a distance r in case 
of Pauli blocked couples Pi (red points, neutron and pro- 
ton couples with same spins), for the case of un-blocked 
proton and neutron couples Pq (black point) and finally 
for neutron-proton couples P„p (blue points) . 

The calculations are referred to a symmetric portion 
of nuclear matter consisting of a sphere containing 2000 
nucleons at a density of about 0.17/m^. The calculations 
include also the iso-vectorial potential energy (T4 ~ 59 
MeV in this example ). The probability distributions 
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FIG. 2. (Color online) Probability distribution Pi to find 
two identical nucleons at a relative distance r. In the same 
figure Po represents the probability evaluated for neutron and 
proton couples with opposite spin. Finally Pnp indicates the 
probability obtained for neutron-proton couples. 



have been estimated by taking in to account only cou- 
ples of nucleons in which at least one of the two nucle- 
ons is localized inside a smaller sphere with a radius of 
4 fm. In this way we can neglect surface effects on this 
quantity. The error bars indicates the uncertainty related 
to the statistics of the simulations. We remark that in 
all the calculations concerning the present work in the 
range of density explored no cluster formation has been 
detected by our numerical procedure which was applied 
at each time step to check for the existence of aggrega- 
tion processes. Only in the lower limit p — O.Tpo and 
for only part of the time, the formation of a big cluster 
with a mass about 98 % the total one and light particles 
have been observed. From the figure it is clearly seen the 
depletion of the Pi distribution at small distances due 
to Pauli constraining in CoMD. For the other distribu- 
tions an enhancement is produced due to the attractive 
effect of the two-body interaction. In particular we note 
that this enhancement is more pronounced for Pnp as 
compared to Pq; this is clearly an effect due to the iso- 
vectorial forces which generate a stronger attraction for 



neutron-proton couples. 

Now from cq.(35) we can obtain an expression for Sy 
as: 



(4^fT2)3/2 
-p(l±/) 



P{r)exp[-j^]dV 



(37) 



47rcr2 



r koexpl 



r 



v{r)dr (38) 



From the above equations we see that even if the overlap 
integral per nucleon is closely related to the density p, 
it admits correction terms related to the spacial correla- 
tion through the integral /. In general, for function v{r) 
localized within a distance of the order of I the value of 
/ decreases with the ratio x — for example, for an 
exponential profile e^T the correction terms go to zero 

like e . In other words, according to what observed 
in the previous section, when the single particles space 
distribution is enough large, an averaging of the spacial 
correlations effect is obtained and the CoMD functional 
tends to the one represented in eq.(lO) which is typical of 
the semi-classical mean-field approximation in transport 
theory. On the contrary for well localized wave-packets 
with Ur of the order of 1-2 /m, / is different from zero 
and reflects the behavior of v{r) at small distances. 

In the case of asymmetric NM we have to generalize 
the above expressions; to this aim it is useful to start 
from the total overlap S^, already introduced in eq.(21) 



— '^1/2,1/2 '^-1/2,-1/2 •-'-1/2,1/2 •^1/2,-1/2 + 

'Tip 

(39) 



qP 1 qP 

*-^l/2,l/2 '-^-1/2,-1/2 



CP 

'-^-1/2,1/2 



CP I C'lp 

"'1/2,-1/2 + 



In the above expression the different terms indicate the 
partial contributions to the total overlap for neutron- 
neutron, proton-proton and neutron-proton couples with 
different combinations of spins. By assuming equal num- 
ber of nucleons with third spin component 1/2 and -1/2 
and applying the generalization of eq.(37) to the different 
subsystems we get: 



o,, 



PnN 



PnN 



2 



(1 



i) 



2 (l-^l,l) + ^(l-^,l) + (pn^ + PpA^)(l+^") 

(40) 

The quantities Is,t indicate the integral containing the 
related spacial correlation functions given in eq.(37). s 
corresponds to third component of the total spin for the 
considered couples (s =0 or l)and t the related isospin 
component (t =0 or 1). Finally /° = /o,o + 

By expressing the partial density and nucleon numbers 
as a function of A and p we finally obtain for the 



7 



following expression: 



-?[8 + (l + /9)'(/o.-i-/i.-i) 
o 

+ (1 - /3)2(/o,i - + 4(1 - (41) 

with an analogous procedure for the other main quanti- 
ties we get: 



/i.-i) 



PA = ^[4+(l + 2/3-2/33)(/o^_i 



(1 + /0) 



[(1 + 6+(/o,-i - + 6-(/o,i - /i,i)] 



(42) 



1 (43) 



(44) 



The above expressions have been obtained by making 
the following approximation — 1-/3^. The pres- 

ence of odd powers in (3 in eq.(41) does not violate the 
charge-symmetry of our system in fact all the above ex- 
pressions are invariant respect to the exchange of neutron 
with protons and a simultaneous change of sign of the 
/3 parameter. The integrals Is,t can have in general an 
intrinsic dependence on /3 and p due to the complex self- 
consistent dynamics, but the determination of the cor- 
relation functions v{r) and of the related integral / is a 
problem which can not be solved in a general way. Some 
special cases are well known. They concern the study 
of the density fluctuations in the hydrodynamical limit 
valid for large distances compared to the mean free-path 
in non-equilibrium cases (see also [IS [HI)- 

However in 

this limit simultaneous spacial correlations are supposed 
to be zero at small distances. At short distances, com- 
parable with the range of the effective interaction, as in 
our case, approximations schemes can be developed ac- 
tually corresponding to a power decomposition in p (see 
the next sections). 

The above relations (eq. (40-43)) suggest a /? depen- 
dence of the primary quantities which define the to- 
tal energy and, in particular, the /? dependence of Sy 
could produce a further coupling between Iso-vectorial 
and Iso-scalar forces. We note that due to the rather 
high strength of the iso-scalar forces as compared with 
the iso-vectorial one, also small changes of the total over- 
lap integral per nucleon as a function of /3 could affect the 
density behavior of the symmetry energy. Nevertheless 
in the present calculations and for the range of asym- 
metry investigated (/3 = — 0.2), the dependence of Sy 
on /3 results to be rather small and comparable with the 
precision of the calculations as shown in Fig. 3. 

To draw definite conclusions on this delicate aspect, 
simulations in a wider range of j3 and with a number 
of particles several times larger than number of particles 
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FIG. 3. Typical result for the total overlap per nucleon Sv as 
a function of the reduced density computed for two different 
values of the asymmetry parameter j3 as shown in the figure. 



used in this work (see the next section) should be per- 
formed. The current computers performances make still 
really difficult this kind of investigations (the cpu time in 
our model calculations depends on a quadratic way from 
the number of particles). 

Therefore concerning the present study we can assume 
that, within the precision of our calculation as shown in 
Fig. 3, the explicit dependence of Sv on /3 is well com- 
pensated by the intrinsic P dependence of the Is,t spatial 
correlation integrals and therefore in the following we will 
refer to the average value of Sv respect to /3 at the dif- 
ferent densities. 



III. 



THE NUCLEAR MATTER SIMULATION 



As mentioned in the introduction, the main goal of 
the present work is the understanding also at a quanti- 
tative level and for a given simple form of the effective 
interaction, of the consequences produced by the corre- 
lations above discussed by taking as a reference point 
the results obtainable in the framework of the Se-MFA 
approach. The study will be performed in a narrow re- 
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gion of densities around po and for a limited range of /? 
values. To this aim in the following we will try to find 
the set of parameters for the effective interaction which 
reproduce in both the cases some of the commonly ac- 
cepted saturation properties of symmetric nuclear mat- 
ter. In particular we will refer to po=0.165 fm^^, bind- 
ing energy E{po) = —16 MeV, compressibility modulus 
for symmetric nuclear mater Knm{po) = 220 MeV and 
for charge/mass asymmetric NM to a symmetry energy 



value 



28 MeV at the saturation density. In the 



case of the Se-MFA, from the functional E reported in 
eqs(10-12), we can find easily the parameter values satis- 
fying the above conditions by solving the system formed 
by the following set of equations for symmetric NM: 



po = 0.165 
E{pa) = -WMeV 
dE 

~ I po 



— In^ =0 



dp 



9a 



d'^E 



Vpo = 220MeV 



(45) 
(46) 

(47) 
(48) 



In concrete cases due to the finite steps with which we 
perform the variation on the parameters the values of 
E{po) and K^Mipa) are obtained within ±0.5%. Within 
the above specified uncertainty the solution gives the fol- 
lowing values for the parameters: Tq — —263 MeVjTa — 
208 MeV and a = 1.25. The value of has been fixed to 
32 MeV which produce a symmetry energy ea{po) = 28.6 
MeV according to eq.(14). 



III.l. NM calculations and CoMD model 

The evaluation of the total energy per nucleon E'~^ re- 
lated to the CoMD calculations requires the solution of 
the many-body problem using the equation of motion 
regulating the wave-packet dynamics. 

At the different densities changing between 0.7-1.2/9o 
with steps equal to O.lpo, the calculations have been 
performed by enclosing a relatively large number of par- 
ticles Ai = 1600 and A2 = 3560 in spherical volumes 
of radii R = tqA^/'^ with tq = ij:^y^^- Particles try- 
ing to escape from the spheres have been re-scattered 
inside through an elastic reflection at the surface. For 
symmetrical configurations, starting from the parame- 
ter values minimizing the functional E in eq. (10-12), we 
have searched for the stationary minimum energy con- 
ditions by applying the cooling-warming procedure cou- 
pled with the constraint related to the Pauli principle 
Ref . . Calculations have been performed for the two 
systems having number of particles equal to Ai and A2 . 
The value of T4 has been fixed to 32 MeV and the cal- 
culations have been performed for different values of the 
stiffness parameter 7. From the minimum energy config- 
urations we have evaluated the related intensive quanti- 
ties Sv{p), a(p),pA{p) and E^in. Corrections due to the 



0.2 




FIG. 4. (Color online) Typical result for the primary quanti- 
ties Svip), a{p),pA{p) as a function of reduced density com- 
puted for P = and 7 = 1. The star symbols indicate the 
results of the study performed on the lighter system with 
mass Ai, the squares represent the results obtained for heav- 
ier system containing A2 particles. Finally the circles indicate 
the obtained corrected values for surface effect according to 
eq.(48). The black solid lines are the final results of a fit pro- 
cedure with a second order polynomial of the density. The line 
plotted with dots in the upper panel represents the density p 
as a function of the reduced density. 



surface effects, which are necessary to estimate the asso- 
ciate bulk values have been evaluated using the following 
relation: 



Qi = Qb + QsA] 



-1/3 



0A5QsA, 



-2/3 



(49) 



Qi indicates the quantity valuated for the system with 
mass Ai (i=l,2). Qi, and Qs are the bulk and surface 
coefficient. Effects related to the curvature are repre- 
sented by the last term of eq.(49). The coefficient 0.45 
has been deduced performing a couple of calculations in 
boxes having the same volume as the considered spheres. 

As an example, for /3 = and 7=1 (form factor ^2) 
in Fig. 4 we show as a function of the reduced density the 
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values of Sv{p), a{p), pa{p) evaluated for the systems 
1 (marked points with star), for the system 2 (marked 
points with square) and the bulk estimated values (dot 
points). Corrections of the same order are evaluated also 
for Ekin (see Sec. II. 2). In the following for simplicity 
we will refer to the bulk quantities without using the 
subscript h. After this first step, we fit the evaluated 
bulk quantities with a polynomial function of the reduced 
density The above quantities show in fact deviations 
from a simple linear behavior as a function of p. This can 
be seen by looking at the red line in the figure (p as a 
function of p/po)- A second order polynomial reproduces 
very well the behavior in the range of explored densities. 
The results of the fit are shown in the figures with lines. 

The obtained functions are then substituted in E'~^ and 
therefore the total binding energy can be now evaluated 
with continuity as a function of p. We can therefore 
search for the parameter values solving the CoMD func- 
tional E'-^ obtained in Sec. II. 2 and satisfying the condi- 
tions expressed in eqs. (45-48). We have to note that the 
numerical solution of this system of coupled equations 
in general can not be obtained with the same precision 
as the one involving the functional E. In most of the 
cases, depending on the stiffness parameter 7 it is possi- 
ble to obtain solutions reproducing the -E(po) and po val- 
ues within 10% while a larger spread is obtained for the 
KNAiipo)- The chosen solutions will be the ones which 
minimize the total relative difference from the reference 
values. Having found the best solution for the functional 
E'-^ , in the sense above specified, with the new set of 
parameter values for Tq, T3 and a we perform another 
series of microscopic NM simulations on the systems of 
Ai,A2 particles. After having included the correction for 
surface effects, we do the polynomial fit. By using the 
new calculated quantities we solve another time the func- 
tional E"-^ and we obtain a new set of parameter values. 
This iterative procedure is continued until the values of 
the parameters differs, in two subsequent steps, by an 
amount less than ±5 %. The method rapidly converges 
after 2-3 iterations. 

After having found the set of parameters values repro- 
ducing (within the above specified precision) the satura- 
tion properties of symmetric NM, for different values of 
the 7 parameter, with others numerical simulations we 
study the system for /3 ^ 0, (/3 ranges from to 0.2 with 
a step 0.05) and we apply the cooling-warming procedure 
until stationary conditions are reached. Through this last 
stage, after the usual corrections for surface effects and 
the polynomial fit as a function of p, the value of E"-^ can 
be computed also for asymmetric NM. This will allows 
us to evaluate the bulk symmetry energy and the related 
density dependence as produced by the model. 



IV. RESULTS AND DISCUSSION 

In this section we illustrate the results obtained from 
the recursive procedure previously described. 




0.7 0.8 0.9 1 1.1 1,2 1.3 
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FIG. 5. (Color online) Final values of Sv{p) for /3 = 0.05 as 
a function of reduced density and for different form factors 
Fk as indicated in the legend. The black lines are the result 
of the fit procedure with a second order polynomial of the 
density. The solid line represents the density p. 

IV. 1. Results on the primary quantities 

As an example in Fig. 5 we show the final values of 
Sy as a function of the reduced density obtained in the 
case of /3 = 0.05 and for different form factors Fk. The 
lines represent the results of the fit with a second order 
polynomial. The red line represents instead the linear 
relation corresponding to the density p. As we can see 
Sy shows deviations from p depending also on the used 
form factor. 

Under the same conditions in Fig. 6 with black line and 
points, we show the value of a as a function of the reduced 
density. The red points represent the values of a obtained 
in the case of T4 = MeV and (3 = using the form 
factor F2. The blue ones represent the values obtained 
for T4 = 59 MeV which corresponds in the case of the 
Se-MFA to a value of esym(po) equal to about 42 MeV. 

As can be seen the value of a decreases with the density 
and increases with T4. For r4 = finite values of a 
are essentially due to correlations imposed by the Pauli 
principle in the system of interacting particles. 
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FIG. 6. (Color online) Final values of the a correlation co- 
efficient for /3 = as a function of reduced density and for 
different form factors Fk as indicated in the legend. The as- 
terisk and full dot symbols (red and blue online) joined with 
continuous lines are associated to T4, values and 59 MeV 
respectively. The other ones are associated to Ta — 32 MeV. 
The solid lines represent fit results with a second order poly- 
nomial of the density. 



IV. 2. The case of symmetric NM 

In Fig. 7, upper panel, we show with a solid line and 
for /3 = the total energy per nucleon as a function of 
the reduced density (eqs.(ll,12)) satisfying the requested 
conditions on NM saturation properties. This result rep- 
resents our reference point. 

The curves with broken lines represent the results of 
our NM simulations with CoMD model at the first step 
of the iteration for the different indicated form factors 
Ffe. In this case the parameter used for the effective in- 
teraction are the same like the ones obtained from the 
minimization of the energy functional related to the Se- 
MFA. As can be seen the curves in this cases are rather 
different. In particular for Fi and F2 the minimum en- 
ergies are shifted to lower values and the compressibility 
modulus also shows large deviations compared to the cho- 
sen reference value. F3 shows instead even a maximum 
at the saturation density. In the lower panel we show the 
corresponding values of Etias (see eq.(33 )) which are 
proportional to a and T4. Together with the density de- 
pendence of the primary quantities above described, this 
term is the main co-responsible of the observed devia- 
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FIG. 7. Upper panel: total energy per nucleon E for symmet- 
ric NM as a function of the reduced density. The solid line 
represents the results obtained through the Se-MFA. With 
discontinue lines we plot the results for CoMD-II calculations 
corresponding to the first step of the iterative procedure (see 
the text). Different discontinue lines represent results related 
to different form factors Fk according to the legend. Bottom 
panel: for the same parameter values the values of Etias (see 
eq.(33)) are plotted as a function of the reduced density. 



tions. The figure shows that the density dependence of 
Ebias strongly depends at a quantitative and qualitative 
level on the used form factors F^. In particular we note 
an increasing slope from positive to negative values with 
the increasing of the degree of stiffness 7 characterizing 
the behavior of the iso-vectorial forces. We can expect 
therefore that the necessary corrections on the parameter 
values of the Iso-scalar effective interaction to reproduce 
the reference properties of the symmetric NM will show 
a dependence on 7. 

In Fig. 8 we show analogues results after that the self 
consistent iteration procedure has been completed as de- 
scribed in Sec. Ill . In the interval 0.85 < 7 < 1.5 the 
saturation density and binding energies are reproduced 
within some percent. The compressibility is instead ob- 
tained within about 20%. In the figure we show some 
example of these solutions with black lines. In general 
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' 7=0.85 K„=271 MeV 
7=1,75 K..= 532 MeV 



F, K,„=219 MeV 
7=0.75 K„„=465 MeV 
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FIG. 8. (Color online) Total Energy per nucleon obtained 
through CoMD-II calculations for symmetric NM as a func- 
tion of the reduced density. Different curves refer to different 
form factors Fk- In the legend the values of Knm at po are 
also indicated. The thin dashed and the thick continues lines 
(red and magenta online) represent an example of limiting 
cases reached when 7 is outside the range 1.5-0.875. 



we note a more pronounced asymmetry of the curves 
around the saturation density with a reduced slope of the 
lower density branch and an increased slope for reduced 
density larger than 1. This behaviors is a direct conse- 
quence of the density dependence of the average overlap 
integral (see Fig. 5). For stiffness parameter values out 
of the indicated interval we observe a fast increasing of 
the compressibility (beyond 400 MeV) and of the corre- 
sponding binding energy at the saturation density. In 
Fig. 8 an example showing this trend is represented by 
colored lines obtained for 7 = 0.75 and 7 = 1.75. All 
these circumstances shows that, in the frame work of 
the present molecular dynamics approach the form 
factor with stiffness parameter changing in the range 
7 ~ 0.5 — 0.85 the due to the correlations discussed is 
not able to reproduce the commonly accepted saturation 
properties of the symmetric nuclear matter. These re- 
sults are consistent with recent findings on the study of 
the 40.48(70+40,48 systems at 25 MeV/nucleon 
concerning the yield balance between incomplete-fusion 
and multi-breakup processes. Finally in Fig. 9 we show 
as a function of the 7 parameter the values of Tq, T3 
and cr as obtained from the iterative procedure described 
in Sec. III.l. Apart from the extremal plotted values 
corresponding to high values of the compressibility ( see 



> 




FIG. 9. The values of the parameters To T3 and a obtained 
through the iterative procedure applied to CoMD-II calcu- 
lations (see the Sec. III.l) are shown as a function of the 
stiffness parameter 7. The values of 7 equal to 1 and 1.5 are 
associated to the functional form F2 and F\ respectively (see 
Sec. II). The others ones are instead associated to F4. The 
bar errors represent the estimated global uncertainty on the 
parameter values related to the numerical procedures. 



Fig. 8), the internal ones correspond to saturation den- 
sity, binding energies and Kj^m values within 15% of 
the value obtained in the case of the Se-MFA. From the 
figure we observe a dependence of the parameter values 
describing the iso-scalar forces on the stiffness parame- 
ter 7 associated to the iso-vectorial interaction. In the 
internal region of the explored interval, even if the de- 
pendence is moderate, maximum change of the order of 
16% and 20% are obtained for the T3 and a parameters 
respectively. However, it is remarkable that the average 
values of the iso-scalar interaction parameters show large 
differences compared to the reference values obtained in 
the case of the Se-MFA (Tq = -263 MeV,T3 = 208 MeV 
and a — 1.25. see Sec.III)which are independent on 7. 



IV. 3. The symmetry energy 

According to what observed in Sec. II. 2 and to 
eq. (14,32), for the range of asymmetries investigated we 
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FIG. 10. Efy^ as a function of /3 is shown for the F2 form 
factor at different densities. The lines represent the results of 
fits with a dependence. 



FIG. 11. (Color online) Panel (a): as a function of the 

density is plotted for the indicated form factors. Panel (b): 
esym in the case of the Se-MFA. 



can express the symmetry energy at different density as: 

Eg^ {p, /3) = E"" {p, P) -E^{p,0) = 
Eg^ (p, /3) - Eg^ {p, 0) + Ek^n if, /3) - Ek^n {p, 0) (50) 
and 



(51) 



The first two terms in eq.(50) represent the change of the 
iso- vectorial interaction while the last two are associated 
with the kinetic energy change. This last variation, after 
the correction for surface effects, according to the con- 
straint is well described by the one related to the Fermi 
motion (see Sec. II. 2). The iso-vectorial interaction con- 
tribution contains the effects related to the discussed cor- 
relations. As an example for different densities we plot 
in Fig. 10 the behavior of E^y^ as a function of (3 for the 
form factor i^2- The dependence well fit the obtained 
values as shown through the lines. 

In Fig.ll. we show in the panel (a) the values of e^,„ 
as a function of the density for different form factors. 
They can be compared with the corresponding values 
in Se-MFA by looking at the panel (b). Even if they 
have similar behavior some remarkable differences are 
obtained in the values of the slope and and curvature 



related to the density dependence. In particular in the 
case of CoMD calculations around the saturation density, 
lower values (of about 1 MeV) are obtained as results of 
the different structure of the iso-vectorial term in eq.(32). 
To investigate in more details on the density dependence 
of the symmetry energy coefficient e^„, as tradition- 
ally done[26l[27|. we perform a Taylor expansion of 
around the saturation density up to the second order: 



-'sym \ 



er„^(p) = e. 



Apq) 



3^ Po 



Po 



Po. 
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PO 



(52) 



L and Ksym are proportional to the slope and curvature 
associated to the density dependence. They determine 
both the corrections on the pressure and on the com- 
pressibility due to the iso-vectorial forces. In Fig. 12, we 
plot for different values of the 7 parameter as indicated 
with numbers in the panel, the values of Ksym versus L 
(black lines) . The lines in the bottom panel represent the 
corresponding values in the case of Se-MFA. 

As we can see the largest changes are observed in Ksym ■ 
These changes reflect the behavior as a function of den- 
sity of the primary quantities S^^ pA and a appearing 
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FIG. 12. (Color online) The values of K sym are plotted versus 
L for the CoMD results (black line joining full sots) and in 
the case of the Se-MFA (blue line joining star symbols). The 
small numbers near the points represent the value of the 7 
parameter. The solid arrows joining the points are plotted 
only to evidence an approximated trend. 



in the expression which determine the behavior of E^^^ 
(see eq.(32)). The density dependence of these quantities 
which contains the correlation effects discussed in Sec. II 
have indeed finite and positive values of the curvature 
around po as can be seen by looking at Fig. 4. This be- 
havior, on turn, derives from the characteristic way in 
which the overlap between the wave-packets (Gaussian 
in our case) changes with the density and therefore is 
strictly linked with the shape of the wave functions used 
to describe the single nucleons. 



V. SUMMARY AND CONCLUDING REMARKS 

Many-body correlations produced in Molecular Dy- 
namics approach based on the CoMD model have been 
discussed and their connection with the used effective in- 



teraction have been analyzed in the case of asymmetric 
nuclear matter simulations. This study has been per- 
formed by comparing the results obtained for the total 
energy functional, and the NM main saturation proper- 
ties with the ones obtainable in the case of a semiclassi- 
cal mean field approximations (Se-MFA). This compari- 
son has been performed by using the same kind of sim- 
ple effective Skyrme interaction. While we can expect 
the two approaches to produce large differences at low 
density due to the cluster formation process, the follow- 
ing study shows that noticeable differences are obtained 
also in a narrow range of densities around the satura- 
tion one where no cluster production has been observed. 
The effective Skyrme interaction include two-body and 
three-body effective iso-scalar interactions plus a two- 
body iso- vectorial one with form factors commonly used 
also in Se-MFA. In CoMD calculations the effects related 
to the spacial correlations produced by the usage of the 
localized wave-packets and the ones related to the multi- 
particles correlations produced through the Pauli prin- 
ciple constraint have been expressed through a second 
order polynomial decomposition of the total energy as a 
function of the density. The obtained results show that, 
contrary to the case of the Se-MFA, the discussed corre- 
lations produce an interdependence between parameters 
describing the iso-scalar forces and the ones related to 
the iso-vectorial interaction. The usage of an iterative 
procedure tuned to obtain in both the cases (Se-MFA 
and in CoMD approaches) very similar saturation den- 
sity, total energy and compressibility for symmetric NM 
(fixed to commonly accepted values) allowed to extract 
a " good" set of parameters used for CoMD calculations. 
The value obtained differs under many aspects from the 
ones obtained from the Se-MFA: in particular the density 
dependence of the used form factors describing the iso- 
vectorial forces can change now in rather more restricted 
range of stiffness values to reproduce the above NM prop- 
erties. Moreover, the values of the coefficients describing 
the iso-scalar interactions are rather different in the two 
cases. For asymmetric NM, in the range of the investi- 
gated asymmetry parameters (/3 = 0. — 0.2), by using the 
same strength for the iso-vectorial interaction, the values 
of the obtained symmetry energies around the satura- 
tion density are only slightly different, but rather large 
differences are instead obtained for the slope L and es- 
pecially for the Ksym curvature parameters. In CoMD 
calculations Ksym assume in fact r higher values. This 
result can be attributed to the shape of the single par- 
ticle wave-packets which determines the behavior of the 
average overlap per nucleon as a function of the density. 
Finally we conclude by observing that even if from a nu- 
merical point of view the obtained results are strictly 
valid for the CoMD model, the performed study shows 
that the observed differences in the parameter values de- 
scribing the effective interaction can have a more wide 
meaning. In fact, they are strictly linked to some general 
properties of the semiclassical wave packets dynamics. 
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